Investigation of molecular penetration depth variation with SMBI fluxes
Zhou Yu-Lin1, Wang Zhan-Hui1, †, , Xu Min1, Wang Qi2, Nie Lin1, Feng Hao3, Sun Wei-Guo4, ‡,
Southwestern Institute of Physics, Chengdu 610041, China
School of Physics and Optoelectronic Engineering, Dalian University of Technology, Dalian 116023, China
School of Physics and Chemistry, Xihua University, Chengdu 610041, China
Institute of Atomic and Molecular Physics, Sichuan University, Chengdu 610041, China

 

† Corresponding author. E-mail: zhwang@swip.ac.cn

‡ Corresponding author. E-mail: swg@mail.xhu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11375053, 11575055, 11405022, and 11405112), the Chinese National Fusion Project for ITER (Grant Nos. 2013GB107001 and 2013GB112005), the International S&T Cooperation Program of China (Grant No. 2015DFA61760), and the Funds of the Youth Innovation Team of Science and Technology in Sichuan Province of China (Grant No. 2014TD0023).

Abstract
Abstract

We study the molecular penetration depth variation with the SMBI fluxes. The molecular transport process and the penetration depth during SMBI with various injection velocities and densities are simulated and compared. It is found that the penetration depth of molecules strongly depends on the radial convective transport of SMBI and it increases with the increase of the injection velocity. The penetration depth does not vary much once the SMBI injection density is larger than a critical value due to the dramatic increase of the dissociation rate on the fueling path. An effective way to improve the SMBI penetration depth has been predicted, which is SMBI with a large radial injection velocity and a lower molecule injection density than the critical density.

1. Introduction

Plasma fueling with higher efficiency and deeper injection is very crucial to meet the fusion power performance requirements for the next generation devices, like ITER. Three main fueling methods, pellet injection (PI),[1] gas puffing (GP),[2] and supersonic molecular beam injection (SMBI),[3] are usually applied for tokamak plasma fueling. The pellet injection is very expensive and requires complex equipment, even though the condensed hydrogen pellets can deeply penetrate into the core plasma with a high injection speed. The fueling efficiency of GP is very low because of its shallow penetration depth. As a new fueling method, SMBI enhances the particle directed group injection speed compared to that of GP (i.e., spreading in all directions with thermal velocity only), and can get a deeper penetration and a higher fueling efficiency than GP. The fuelling efficiency of SMBI is in the range 30%–60%, which is about three times as high as that of GP. In addition, its cost is much lower than that of PI. SMBI has been widely used in tokamaks such as HL-2A,[4] EAST,[5] KSTAR,[6] and JT-60U.[7,8] Besides being a fuelling method, SMBI is also a good tool for the L–H transition and confinement improvements,[9] edge localized mode (ELM) mitigation,[10] nonlocal heat transport,[11] and impurity transport[12] studies. Until now, the mechanism causing a deeper penetration and a higher fueling efficiency of SMBI has not been well understood yet. Thus, researches on the SMBI fueling efficiency depending on different fueling fluxes have been promoted in both experiments and large scale simulations.

The SMBI fueling efficiency mainly depends on the penetration depth, which is affected by the different fueling fluxes. The studies of the different fueling fluxes on the penetration depth are very important to understand the neutral particle transport dynamics and improve the fueling efficiency of SMBI. Many experiments have been carried out to study the mechanism of the high fuelling efficiency of SMBI. For instance, the fuelling efficiency was compared between GP and SMBI on Tore Supra (TS).[13]

Owing to the technical problems and the complicated plasma environment, lots of crucial physical processes cannot be directly measured in the experiment with high temporal and spatial resolutions, such as the convection reversal due to the turbulence, the penetration front of neutrals, the 2D/3D transport dynamics of particles and their atomic collisional interactions. The large scale simulation is a more efficient and economical way to study these problems. For example, a fluid code based on the ion temperature gradient (ITG) mode and the trapped electron mode (TEM) has been developed recently to study the spatial-temporal dynamics of the anomalous particle convection reversal.[14] A gyro-fluid code has been developed and used to study the interaction between two-dimensional vortex flows and micro-turbulence,[15] and the characteristics of the ion temperature gradient (ITG) instability in the presence of a magnetic island.[16] Thus, it is also necessary and helpful to study and understand the fueling mechanism with the reduced transport model or large scale simulations. Since the kinetic particle collision simulations with the Monte–Carlo method demand a great amount of computer time and resources, it is more efficient to simulate the neutral transport and interaction with plasma with the fluid method during SMBI (i.e., high molecule injection density and high collision reaction rates with plasmas).

Recently, the BOUT++ code has included a new module, named trans-neut,[17] to study the time-dependent 3D transport dynamics during SMBI in real tokamak geometry. With this module, the time-dependent transport of both plasma and neutrals was simulated during SMBI, including the front propagation of neutrals and the evolution of the edge plasma and neutral profiles.[17] It revealed the transport dynamic in both radial and poloidal directions during SMBI. The numerical initial profiles of the plasma density and temperatures were selected to be the same as the experiment profiles right before SMBI at 830 ms of the shot # 22535 in the HL-2A tokamak.[18] The total effective injected particles (the ones finally being ionized to become plasma) in simulations were the same as those in the experiment. It has been found that the variations of the plasma profiles during SMBI consist qualitatively well with each other between the simulation and the experiment.[18] With this module, it has also been found that the radial convection and larger inflowing flux compared to GP lead to the deeper penetration depth and the higher fueling efficiency of SMBI.[19] There were some other approximate approaches developed by Tokar et al. in Refs. [20]–[22] which reduce the three-dimensional transport equations to a set of one-dimensional ones. There were also some other codes to study the neutral transport and interactions with plasma during GP, such as UEDGE,[23,24] TOKES,[25] DEGAS,[26] and SOLPS[27] but just few modeling programs and codes were able to simulate SMBI.

This paper is organized as follows. The physical model is simply presented in Section 2. Section 3 gives the simulation results of transport dynamics and neutrals penetration during fueling with different injection velocities and injection densities. Finally, these results are summarized in Section 4.

2. Physical model
2.1. Transport equations

The plasma density, heat, and momentum transport equations together with the neutral density and momentum transport equations are included in a seven-field fluid model. Some basic physical effects are considered during the fueling process, including the perpendicular plasma density diffusion, heat diffusion, parallel plasma density convection and heat conduction, energy interchange between ion and electron, parallel ion viscosity, atom diffusion, and molecule radial density convection. From Ref. [17], the transport equations are re-written as

Equations (1)–(4) describe the plasma density, heat, and momentum transports, while equation (4) describes the atom density transport, and equations (5) and (6) denote the molecular density and momentum transports. The plasma quasi-neutral condition Ni = Ne is applied in the physical model. In these equations, , , and are the perpendicular classical diffusion coefficients of ion density, ion temperature, and electron temperature, respectively. and are the perpendicular (to the magnetic field) and parallel atom diffusion coefficients. and are the parallel classical thermal conductivity coefficients. is the atom ionization source, and is the ion and electron recombination sink of plasma density. νI, νrec, νdiss, and νCX are the atom ionization rate, plasma recombination rate, molecule dissociation rate, and ion–atom charge exchange rate, respectively. Wrec is the fraction of recombination energy re-absorbed by an electron during recombination via some processes (i.e., three-body recombination); WI and Wdiss are the electron energy lost per ionization and dissociation; and Wbind is the binding energy between the two hydrogen atoms of a molecule. is the parallel ion viscosity. P and Pm are the total plasma pressure and molecule pressure where the molecule temperature Tm is the room temperature (i.e., 300 K). ∇x is the radial component of the gradient. The specific details can be found in the reference.[17]

2.2. Boundary conditions

In the HL-2A tokamak, the poloidal cross-section of the tokamak geometry with X-point was separated into three principle regions in Refs. [28] and [29]: the plasma core and edge region, the scrape-off-layer (SOL) region, and the private flux and divertor plate region.[30] In the simulations, the fueling source is injected to the low field side of the HL-2A tokamak.

In tokamaks, plasma fueling like GP and SMBI is to inject molecules with a constant inflowing particle flux localized in both poloidal and toroidal directions at the edge. It is thus more realistic to simulate this fueling process by applying the constant flux boundary conditions locally at the edge than giving a simple source term in the equations. Therefore, the local constant flux boundary conditions of molecular injection density and velocity are described as follows:

In Eqs. (8) and (9), the whole fueling range of the toroidal simulation domain is localized within a poloidal range (θ0θθ1), while the fueling source center is at θ1/2 = (θ0 +θ1)/2. wθ is the width of fueling at the outermost boundary flux surface and a = 40 cm is the minor radius.

Besides, the radial (x) flux-driven boundary conditions at the innermost and the outermost boundary flux surfaces, the sheath boundary conditions at the divertor plates, and the particle recycling boundary conditions on both wall and divertor plates have also been included in the simulations. More details about the boundary conditions can be found in the reference.[17]

3. Results

To understand the mechanism of high fueling efficiency and deep penetration, the SMBI pulses under different injection fluxes are simulated and compared below by varying the injection velocity and the injection density. First of all, it is simulated and compared in Subsection 3.1 with the same molecule injection density Nm = 0.5N0 (N0 = 1019 m−3) but different injection velocities varying from Vm = 160 m/s to Vm = 3000 m/s. Then it is simulated and compared in Subsection 3.2 with the same SMBI velocity (Vm = 800 m/s) but different molecule injection densities varying from Nm = 0.1N0 to Nm = 5.0N0. All the SMBI pulses are initially injected at the same steady state treated as t = 0 ms, as shown in Fig. 1, and the durations of the SMBI pulses are the same, that is t = 1ms. What is more, the SMBI pulses with the same injection flux (Γ = Nm Vm = 800N0 m/s) but different injection velocities and densities are simulated and compared in Subsection 3.3. To achieve the same injection flux, the injection velocity is varying from Vm = 200 m/s to Vm = 3200 m/s, while the injection density is varying from Nm = 4.0N0 to Nm = 0.25N0. The simulation results are given as follows.

Fig. 1. The evolving quantities of (a) Ti, (b) Te, and (c) Ni at steady state in the radial direction before SMBI.
3.1. Penetration depths varying with SMBI velocities

First, we study the dynamics of neutrals and plasma transport in the radial direction during SMBI with different injection velocities but the same injection density. Figure 2(a) shows the evolution of the radial molecule propagation with the molecule injection velocity varying from Vm = 160 m/s to Vm = 3000 m/s, while the injection densities are the same, i.e., Nm = 0.5N0. It is obvious that the depth of the molecule propagation front is increased through enlarging the molecule injection velocity. When the molecular injection velocity is very small such as Vm = 160 m/s, most of the injected molecules cannot even penetrate into the separatrix. The inward propagation of molecules can be conspicuously observed with the increase of the radial molecular injection velocity. As shown by the red curve in Fig. 2(a), the molecule propagation front can penetrate deeper inside (ψ = 0.8) when Vm = 3000 m/s.

Fig. 2. The radial profiles of (a) Nm, (b) Ti, (c) Te, and (d) Ni along the fueling path at t = 0.3 ms after SMBI with different injection velocities but the same injection density Nm = 0.5N0.

Figures 2(b) and 2(c) are the radial profiles of plasma temperatures. It is clear that the plasma along the fueling path will lose energy during fueling due to the collision and reactions with neutrals, such as ion–atom charge exchange, dissociation, and ionization. The plasma will be cooled wherever the molecule front propagates. The deeper the molecule front propagates with a larger injection velocity, the lower the plasma temperature will be. In Fig. 2(d), the peak of Ni due to the largest atom ionization rate moves inwards with the increase of the injection velocity. The dissociation rate will also be enhanced with the increase of the injection velocity and the penetration depth. Most molecules are dissociated and then ionized to be plasma in the molecular propagation front region where the dissociation rate is the highest as shown in Fig. 3. The molecule dissociation rate is gradually increasing and also moving inwards with the increase of the injection velocity. It is consistent with the inward propagation of the molecule front. The dissociation rate at the front of the injection beam is crucial to the plasma transport process which will balance with the SMBI injection rate and prevent the molecule to penetrate deeper.

Figure 4 describes the time evolution of the parallel plasma density at ψ = 1.05 in the SOL region during SMBI with the same injection density Nm = 0.5N0 and different velocities. The asymmetry of the plasma densities at θ = 0 and θ = 2π is due to the parallel sheath boundary condition on the divertor plates. The plasma density is gradually spreading on the poloidal direction with the drive of the parallel convection. As the parallel velocities (Fig. 5) are varying similar to each other at the same time, the parallel convections are about the same with increasing injection velocities, and the plasma density spread out of the fueling path with a similar rate in the poloidal direction. Due to the increase of the SMBI injection flux, the plasma density increases in the fueling path as shown in Fig. 4. The plasma density grows very slowly at the edge. The reasons are given as follows. 1) The main region of dissociation (Fig. 3) moves deeper from the edge region due to the increasing penetration depth, so most of the injected molecules are dissociated and then ionized at deep inside the separatrix instead of the edge region. 2) The molecular densities (Fig. 2(a)) are similar in the edge region. It causes the molecular dissociation rate, which is proportional to the molecular density, to be about the same at the edge. Thus, the increasing rates of the plasma density are low and similar to each other in the edge region, which are smaller than those in the cases of SMBI with the same velocity but different densities.

The plasma accumulated in the fueling path will collide with neutrals and reduce its penetration speed, thus it is crucial to prevent the plasma density from dramatically growing in the fueling path and then we can get a better penetration depth.

Fig. 3. The time and spatial evolutions of neutrals dissociation under the same density Nm = 0.5N0 but different velocities of (a) Vm = 160 m/s, (b) Vm = 400 m/s, (c) Vm = 800 m/s, (d) Vm = 1600 m/s, (e) Vm = 3000 m/s.
Fig. 4. Time evolution of parallel plasma density at ψ = 1.05 under the same density Nm = 0.5N0 but different velocities of (a) Vm = 160 m/s, (b) Vm = 400 m/s, (c) Vm = 800 m/s, (d) Vm = 1600 m/s, (e) Vm = 3000 m/s.
Fig. 5. Time evolution of parallel ion velocity at ψ = 1.05 under the same density Nm = 0.5N0 but different velocities of (a) Vm = 160 m/s, (b) Vm = 400 m/s, (c) Vm = 800 m/s, (d) Vm = 1600 m/s, (e) Vm = 3000 m/s.
3.2. Penetration depths varying with SMBI densities

There is another way to study the dependence of the penetration depth on the SMBI flux by varying the molecule injection density. Figure 6(a) shows the evolution of the molecular propagation front with the same injection velocity but different injection densities in the radial direction at 0.3 ms after turning on the SMBI. It is apparent that the molecule density decreases almost linearly from the edge towards the core for a relatively small injection density Nm = 0.1N0 or a small injection flux (the black curve in Fig. 6(a)), and the penetration is much deeper than that with a very large injection density Nm = 5.0N0 or a large injection flux (the red curve in Fig. 6(a)). It is found that the penetration depth of SMBI cannot be effectively improved once the molecular injection density is larger than a critical value of Nm = 0.5N0. Moreover, there will be a very sharp gradient region at the propagation front of the molecular beam due to the dramatic increase of the local plasma density and the dissociation rate. Once the molecular injection density is larger than the critical value, the injection beam front of SMBI is always around and moves inwards inconspicuously with the increase of the molecular injection density. The molecule inward propagation is terminated and the fueling path is blocked due to the dramatic increase of the dissociation rate. Once the molecular injection density is lower than the critical value, the plasma density continually spreads out in the parallel direction (i.e., parallel density flux) and does not accumulate in the fueling path to terminate the inward propagation. However, the local plasma density will be accumulated once the molecular injection rate and the dissociation rate are much larger than the parallel plasma density spreading rate. The critical molecular injection density to get a larger penetration depth depends on whether the plasma density accumulates and the dissociation rate dramatically increases on the fueling path to terminate the molecule inward front propagation. Once the local plasma density increasing rate (i.e., atom ionization rate) is smaller than the local decreasing rate (i.e., parallel and radial spreading rate), the plasma density will not accumulate dramatically, otherwise, it will accumulate dramatically and then restrain the molecular penetration depth. Figures 6(b) and 6(c) present the radial profiles of ion and electron temperatures by varying. The plasma temperature decreases where the molecules penetrate due to the strong cooling effects, and the temperature gradient at the beam front is larger with a larger molecular injection density. In Fig. 6(d), the peak increases dramatically with the increase of the molecular injection density, especially for Nm = 5.0N0, which leads to strong molecular dissociation and also the formation of the sharp gradient region at the molecular beam propagating front.

Fig. 6. The radial profiles of (a) Nm, (b) Ti, (c) Te, and (d) Ni at 0.3 ms of SMBI with the same injection velocity (Vm = 800 m/s) but different injection densities.
Fig. 7. The time and spatial evolutions of neutrals dissociation under the same injection speed Vm = 800 m/s but different densities of (a) Nm = 0.1N0, (b) Nm = 0.25N0, (c) Nm = 0.5N0, (d) Nm = 1.0N0, (e) Nm = 5.0N0.

In Fig. 7(a), the molecular dissociation region with a small injection density (Nm = 0.1N0) is deep inside the separatrix. In Figs. 7(b)(7e), most of the injected molecules are dissociated at the front of the injection beam where the dissociation rate is the highest. Due to the fast and great increasing plasma densities in the fueling path (Fig. 6(d)), the largest dissociation rate increases greatly from about 2.5 × 104N0 s−1 to 2.0 × 107N0 s−1 with increasing injection densities. When the total dissociation rate becomes larger enough than the SMBI injection rate, the penetration depth of SMBI will be decreased. In Fig. 8, the plasma density in the fueling path dramatically rises at the edge (i.e.,) when Nm increases, especially for Nm = 5.0N0. The parallel plasma velocities (Fig. 9) are similar at the same time with different injection densities due to the parallel pressure gradients being about the same. With the similar parallel convection, it provides the similar parallel spreading density flux to reduce the local plasma density in the fueling path for different injection densities. Once the local plasma density increasing rate (i.e., atom ionization rate) is larger than the parallel density spreading rate (i.e., parallel spreading density flux), it will lead to the accumulation of the local plasma density in the fueling path. The plasma accumulated in the fueling path will restrict the penetration of neutrals due to the intense collision reaction with neutrals.

Fig. 8. Time evolution of parallel plasma density at ψ = 1.05 under the same injection speed Vm = 800 m/s but different densities of (a) Nm = 0.1N0, (b) Nm = 0.25N0, (c) Nm = 0.5N0, (d) Nm = 1.0N0, (e) Nm = 5.0N0.
Fig. 9. Time evolution of parallel plasma velocity at ψ = 1.05 under the same injection speed Vm = 800 m/s but different densities of (a) Nm = 0.1N0, (b) Nm = 0.25N0, (c) Nm = 0.5N0, (d) Nm = 1.0N0, (e) Nm = 5.0N0.
3.3. The penetration depth for the same injection flux

Both neutral and plasma transport processes during SMBI under the same injection flux are simulated and compared by varying the injection density and the injection velocity. Figure 10(a) shows the radial profiles of molecular density along the fueling path at 0.25 ms after turning on the SMBI with the same injection flux (Γ = Nm Vm = 800N0 m/s). The penetration depth of the molecular propagation front is obviously increasing with the increase of the injection velocity and the decrease of the injection density for the same injection flux. The propagation front of the injection beam with the largest injection velocity Vm = 3200 m/s (the red curve in Fig. 10(a)) can reach about ψ=0.8, while it just reaches the edge region (ψ = 0.96) with the smallest injection velocity Vm = 200 m/s (the black curve in Fig. 10(a)). With the same injection flux, the radial initial injection velocity dominates the radial convective transport of molecules during SMBI, which leads to a deep penetration depth. Figure 10(b) presents the radial profiles of the molecular dissociation rate along the fueling path at 0.25 ms after turning on the SMBI with the same injection flux. The main molecular dissociation regions are all distributed at the injection beam fronts where the dissociation rates are the highest. When the injection velocity increases, the main dissociation region moves inwards consistent with the inward propagation of the molecule front. Meanwhile, the maximum of the molecular dissociation rate decreases and the dissociation region increases to keep the total molecular dissociation rate balance with the molecular injection rate. As shown in Fig. 10(b), the maximum dissociation rate for a high density Nm = 4.0N0 and a low velocity Vm = 200 m/s(the black curve) locating near the edge is about three times larger than that for a low density Nm = 0.25N0 and a high velocity Vm = 3200 m/s (the red curve), which leads to a shallow penetration depth. The peak of the plasma density is larger and distributed near the edge with a high injection density Nm = 4.0N0 and a low velocity Vm = 200 m/s in Fig. 10(c). Figures 10(d) and 10(e) show the radial profiles of plasma temperature along the fueling path at 0.25 ms after turning on the SMBI with the same injection flux. Despite the fact that the fueling fluxes are the same, the radial cooling areas of plasma temperatures are enlarged with a larger injection velocity due to the deeper penetration. It is found that the penetration depth mainly depends on the injection velocity even for the same injection flux of SMBI, which is the same as the rule for different injection fluxes by varying the injection velocity only. It is because the radial convection transport dominates the molecule penetration rather than the thermal diffusion transport due to the low injection temperature (the room temperature).

Fig. 10. The radial profiles of (a) Nm, (b) (molecular dissociation rate), (c) Ni, (d) Ti, and (e) Te along the fueling path at t = 0.25 ms after SMBI with the same injection flux (Γ = Nm Vm = 800N0 m/s) but different injection densities and velocities.
4. Summary

In this paper, the molecular transport process during plasma fueling of SMBI with different injection velocities and densities has been simulated and compared in realistic divertor geometry of the HL-2A tokamak. It is found that the penetration depth of molecules strongly depends on the radial convective transport of SMBI and increases with the increase of the injection velocity. The principal results are summarized as follows.

The simulation results reveal that the radial directed convection is directly related to the penetration depth of SMBI, while SMBI with the injection density larger than the critical value hardly get a deeper penetration depth. An effective way to improve the SMBI penetration depth has been predicted, which is SMBI with a large radial injection velocity and a lower molecule injection density than the critical density. Some following works can be done to study the SMBI: 1) validate the simulation results with experiments by improving the injection velocity to study the variations of the penetration depth; 2) find a way to increase the penetration depth in the H-mode plasma by doing simulations; 3) improve the simulation model by coupling with a turbulence transport equation to study the turbulence influence on plasma and neutral transport in the edge region during SMBI.

Reference
1Baylor L RJernigan T CCombs S KHoulberg W AOwen L WRasmussen D AMaruyama SParks P B 2000 Phys. Plasmas 7 1878
2Sajjad SGao XLing BBhatti S HAng T 2009 Phys. Lett. 373 1133
3Yao L HZhao D WFeng B BChen C Y 2010 Plasma Sci. Technol. 12 529
4Xiao W WZou X LDing X TDong J QYao L HSong S DLiu Z TGao Y DFeng B BSong X MYang Q WYan L WLiu YDuan X RPan C HLiu Y 2010 Rev. Sci. Instrum. 81 013506
5Zheng X WLi J GHu J SLi J HDing RCao BWu J H 2013 Plasma Phys. Control. Fusion 55 115010
6Xiao W WDiamond P HKim W C 2014 Nucl. Fusion 54 023003
7Fasoli AGormenzano CBerk H L 2007 Nucl. Fusion 47 S264
8Kwon MOh Y KYang H L 2011 Nucl. Fusion 51 094006
9Xiao W WZou X LDing X T 2010 Phys. Rev. Lett. 104 215001
10Xiao W WDiamond P HZou X L 2012 Nucl. Fusion 52 114027
11Sun H JDing X TYao L H 2010 Plasma Phys. Control Fusion 52 045003
12Song S DGao J MFu B Z 2013 Chin. Phys. 22 125201
13Pégourié BTsitrone EDejarnac R 2003 J. Nucl. Mater. 313�?16 539
14Sun T TChen S YWang Z H 2015 Chin. Phys. Lett. 32 35201
15Wang Z XLi Q JDong J QKishimoto Y 2009 Phys. Rev. Lett. 103 015004
16Wang Z XLi Q JKishimoto YDong J Q 2009 Phys. Plasmas 16 060703
17Wang Z HXu X QXia T YRognlien T D 2014 Nucl. Fusion 54 043019
18Wang Z HXu X QXu M201425th IAEA Fusion Energy Conference2014 TH/P7-30
19Zhou Y LWang Z HXu X QLi H DFeng HSun W G 2015 Phys. Plasmas 22 012503
20Tokar M Z 2000 Phys. Plasma 7 2432
21Tokar M Z 2008 AIP Conf. Proc. 1061 94
22Tokar M Z 2012 AIP Conf. Proc. 1504 1210
23Rognlien T DBraams B JKnoll D A 1996 Contrib. Plasma Phys. 36 105
24Rognlien T DRyutov D DMattor NPorter G D 1999 Phys. Plasmas 6 1851
25Landman IJaneschitz G J 2007 Nucl. Mater. 363�?65 1061
26Heifetz DPost DPetravic MWeisheit JBateman G 1982 J. Comp. Phys. 46 309
27Schneider RBonnin XBorrass K 2006 Contrib. Plasma Phys. 46 3ߝ191
28Xu X QUmansky M VDudson BSnyder P B 2008 Commun. Comput. Phys. 4 949
29Umansky M VXu X QDudson BLoDestro L LMyra J R 2009 Comput. Phys. Commun. 180 887
30Sang C FDai S YSun J Z 2014 Chin. Phys. 23 115201